library(dplyr)
library(ggplot2)
library(multiplot)
setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201120")
file=read.csv("more.than.2AShM.in.DC.add.BayesFactor.beta.DC.csv",head=T)
file=file[file$BF_in_DC>10,]

normaldata=select(file,unitID,DC.con.beta0,DC.con.beta0.upper,DC.con.beta0.lower,DC.con.FDR,avsnp150)
normaldata$group=ifelse(normaldata$DC.con.FDR<0.1,"sig","nosig")
normaldata.nosig=normaldata[normaldata$group=="nosig",]
normaldata.sig=normaldata[normaldata$group=="sig",]
normaldata.nosig=arrange(normaldata.nosig,DC.con.beta0)
normaldata.sig=arrange(normaldata.sig,DC.con.beta0)
normaldata=rbind(normaldata.sig[normaldata.sig$DC.con.beta0<0,],normaldata.nosig,normaldata.sig[normaldata.sig$DC.con.beta0>0,])

#normaldata=arrange(normaldata,group,DC.con.beta0)
normaldata$num=1:dim(normaldata)[1]
normaldata$na=1:200
orders=select(normaldata,unitID,num)

p1=ggplot(normaldata,aes(y=DC.con.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.con.beta0.lower,ymax=DC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")


tumordata=select(file,unitID,DC.case.beta0,DC.case.beta0.upper,DC.case.beta0.lower,DC.case.FDR,avsnp150)
tumordata$group=ifelse(tumordata$DC.case.FDR<0.1,"sig","nosig")
tumordata=merge(tumordata,orders,by="unitID")
tumordata=arrange(tumordata,num)
tumordata$na=1:200

tumordata1=tumordata[1:98,]
tumordata1$length.beta0=abs(tumordata1$DC.case.beta0.upper-tumordata1$DC.case.beta0.lower)
tumordata1.sig=tumordata1[tumordata1$DC.case.FDR<0.1,]
tumordata1.nosig=tumordata1[tumordata1$DC.case.FDR>0.1,]
tumordata1.nosig=arrange(tumordata1.nosig,DC.case.beta0,length.beta0)
tumordata1.nosig$num=1:dim(tumordata1.nosig)[1]
tumordata1.sig=arrange(tumordata1.sig,DC.case.beta0,length.beta0)
tumordata1.sig$num=50:98
tumordata1=rbind(tumordata1.nosig,tumordata1.sig)
tumordata1=tumordata1[,-10]

tumordata2=tumordata[169:200,]
tumordata2$length.beta0=abs(tumordata2$DC.case.beta0.upper-tumordata2$DC.case.beta0.lower)
tumordata2.sig=tumordata2[tumordata2$DC.case.FDR<0.1,]
tumordata2.nosig=tumordata2[tumordata2$DC.case.FDR>0.1,]
tumordata2.sig=arrange(tumordata2.sig,length.beta0,DC.case.beta0)
tumordata2.sig$num=(169+11):169
tumordata2.nosig=arrange(tumordata2.nosig,DC.case.beta0,length.beta0)
tumordata2.nosig$num=181:200
tumordata2=rbind(tumordata2.sig,tumordata2.nosig)
tumordata2=tumordata2[,-10]

tumordata=rbind(tumordata1,tumordata[99:168,],tumordata2)
orders=select(tumordata,unitID,num)
tumordata=arrange(tumordata,num)
tumordata$na=1:200

p2=ggplot(tumordata,aes(y=DC.case.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.case.beta0.lower,ymax=DC.case.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")

normaldata=select(file,unitID,DC.con.beta0,DC.con.beta0.upper,DC.con.beta0.lower,DC.con.FDR,avsnp150)
normaldata$group=ifelse(normaldata$DC.con.FDR<0.1,"sig","nosig")
normaldata=merge(normaldata,orders,by="unitID")
normaldata=arrange(normaldata,num)
normaldata$na=1:200
p1=ggplot(normaldata,aes(y=DC.con.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.con.beta0.lower,ymax=DC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")


layout <- matrix(c(1,2), nrow = 1)
multiplot(plotlist=list(p1,p2),layout=layout)

write.table(normaldata,"./normaldata.txt",quote = F,row.names = F,sep="\t")
write.table(tumordata,"./tumordata.txt",quote = F,row.names = F,sep="\t")

normaldata=read.table("./normaldata.txt",head=T,sep="\t")
tumordata=read.table("./tumordata.txt",head=T,sep="\t")
normaldata$na=""
tumordata$na=""
p3=ggplot(normaldata,aes(y=DC.con.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.con.beta0.lower,ymax=DC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")

p4=ggplot(tumordata,aes(y=DC.case.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.case.beta0.lower,ymax=DC.case.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
multiplot(plotlist=list(p3,p4),layout=layout)


order=data.frame(unitID=normaldata$unitID,num=normaldata$num)
CC=read.table("CC.beta0.txt",head=T,sep="\t")
HC=read.table("HC.beta0.txt",header = T,sep="\t")

CC=CC[CC$unitID %in% normaldata$unitID,]
HC=HC[HC$unitID %in% normaldata$unitID,]

CC=merge(CC,order,by="unitID")
HC=merge(HC,order,by="unitID")
p5=ggplot(CC,aes(y=DC.con.beta0,x=num))+scale_x_continuous(breaks = CC$num,labels = CC$na,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=DC.con.beta0.lower,ymax=DC.con.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
